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1.  Technical  Objectives 


The  objective  of  this  three-year  research  project  is  to  qualitatively  and  quantitatively  map  the 
relationship  between  GB  characters  (structure)  and  material  cohesive  strengths  with  fatigue  crack 
initiation  and  SCG  mechanisms  (initiation  and  growth  versus  resistance). 

2.  Final  summary 

Here,  we  present  a  hierarchical  form  of  multiscale  modeling  to  understand  fundamental  issues 
of  deformation  in  Ti.  This  report  presents  our  effort  in  the  following  areas, 

1.  Effects  of  oxygen  on  prismatic  faults  in  a-Ti:  A  combined  quantum  mechanics/molecular 
mechanics  study 

2.  Nano-indentation  and  slip  transfer  (critical  in  understanding  crack  initiation) 

3.  An  extended-finite  element  framework  (XFEM)  to  study  SCG  mechanisms 

4.  Atomistic  methods  to  develop  a  grain  and  twin  boundaries  database  in  a-Ti  to  quantify  and 
develop 

a.  a  crack  initiation  model, 

b.  the  necessary  cohesive  zone  models  for  interfaces  and  single  crystals,  and 

c.  a  crystal  plasticity  model 

Fundamental  to  the  development  of  a  hierarchical  atomistic-continuum  methodology  for  SCG 
is  the  development  of  appropriate  crack  initiation  models,  cohesive  zone  properties,  crystal 
plasticity  models,  and  the  microstructure  based  finite  element  framework,  i.e.,  in  this  case,  we  are 
using  an  XFEM.  Overall,  we  have  generalized  XFEM  to  include  implicit  representation  of  grains 
and  grain  boundaries.  The  new  formulation  can  incorporate  microstructure  features  directly  based 
on  EBSD  data  using  level  set  formulations.  Furthermore,  an  atomistic  database  of  interactions  of 
dislocation-impurities  was  developed.  Using  the  developed  database  an  appropriate  crystal 
plasticity  model  was  developed  to  investigate  SCG  in  Ti  alloys. 


Page  |  1 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


3.  Technical  details 

3.1  Effects  of  oxygen  on  prismatic  faults  in  a-Ti:  A  combined  quantum  mechanics/molecular 
mechanics  study 

Titanium  (Ti)  and  its  alloys  have  been  traditionally  used  for  structural  applications  in 
automotive,  aerospace,  and  biomedical  applications  due  to  their  high  strength-to-weight  ratios  [1], 
Commercially  pure  Ti-alloys  are  particularly  attractive  due  to  their  excellent  corrosion  resistance, 
light  weight,  and  formability.  Crucial  to  the  mechanical  properties  of  these  alloys  is  the  presence 
of  oxygen  impurities.  Specifically,  it  has  been  shown  that  the  oxygen  impurities  can  either  harden 
or  soften  the  materials  as  a  result  of  the  interaction  between  the  impurities  and  lattice  defects, 
particularly  dislocations  and  twins  [2 — 4].  In  high  strength  grade  Ti-alloys  the  addition  of  oxygen 
improves  the  corrosion  and  wear  resistance  of  a-Ti  and  its  alloys  [1,5].  However,  at  elevated 
temperatures  an  oxide  layer  at  the  Ti  surface  is  formed  and  the  oxygen  atoms  can  diffuse  rapidly 
into  the  base  metal  [4]  instigating  hardening  and  embrittlement  of  the  base  material  [3].  In 
particular,  a  dislocation  core  could  act  as  a  fast  path  for  diffusing  oxygen  atoms  whose  mobility 
can  be  orders  of  magnitude  higher  than  those  in  bulk  diffusion,  and  this  phenomenon  is  often 
referred  to  as  “pipe  diffusion”.  By  creating  a  short-circuit  pathway  the  pipe  diffusion  can  affect 
many  kinetic  processes  in  the  bulk  materials,  including  creep,  [6]  dynamic-strain  aging, 
crystallization  [7],  and  other  mechanical  properties  of  a-Ti.  Consequently,  understanding  the 
atomistic  mechanisms  behind  oxygen-metal  interaction  is  crucial  for  developing  high-performance 
Ti  alloys  [2,8,9]  with  applications  in  nuclear  power  plants,  wind  turbines,  and  other  large-scale 
industrial  infrastructure. 

In  this  report  we  examine  oxygen  diffusion  in  bulk  a-Ti  and  pipe  diffusion  along  a  prismatic 
edge  dislocation  core.  In  addition,  we  evaluate  the  dislocation  mobility  in  the  presence  and  absence 
of  oxygen  impurities.  We  have  employed  a  coupled  quantum  mechanics  and  molecular  mechanics 
(QM/MM)  method  that  combines  the  accuracy  of  QM  calculations  with  the  simplicity  and 
efficiency  of  MM  simulations.  In  conjunction  with  a  climb-based  nudged  elastic  band  (C-NEB) 
method  [10]  we  can  determine  the  oxygen  diffusion  barriers  accurately.  More  specifically,  the 
QM/MM  model  is  partitioned  into  two  spatial  domains:  a  QM  region  consisting  of  the  dislocation 
core  and  the  oxygen  impurity  and  an  MM  region  for  the  rest  of  the  system.  In  the  QM  region  the 
energy  and  force  are  calculated  based  on  the  density  functional  theory  (DFT)  [11,12]  as 
implemented  in  the  Vienna  Ab-initio  Simulation  Package  (VASP)  [13]  while  the  same  quantities 

Page | 2 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


in  the  MM  region  are  computed  by  empirical  atomistic  simulations  [14]  based  on  the  embedded 
atom  method  (EAM)  due  to  Zope  and  Mishin  [15].  The  coupling  between  the  QM  and  MM  regions 
is  accomplished  by  the  constrained  DFT  [16].  The  post  processing  of  the  simulation  data  is 
performed  using  OVITO  [17].  Complementary  to  the  atomistic  simulations,  the  Semi-discrete 
Variational  Peierls-Nabarro  model  (SVPN)  [18]  was  also  used  to  analyze  the  dislocation  core 
structure  and  mobility.  The  chemical  interaction  between  the  O  impurity  and  the  dislocation  core 
is  captured  by  the  so-called  generalized  stacking  fault  energy  (GSFE)  surface  which  was 
determined  from  DFT-VASP  calculations.  By  taking  the  chemical  interaction  into  consideration 
the  SVPN  model  can  predict  the  dislocation  core  structure  and  mobility  in  the  presence  and 
absence  of  the  O  impurity  and  thus  reveal  the  effect  of  O  on  the  softening/hardening  behavior  in 
a-Ti. 

An  edge  dislocation  was  created,  as  described  by  Osetsky  and  Bacon  [19],  initially  in 
LAMMPS  [20],  The  edge  dislocation  structure  obtained  using  the  EAM  potential  was  further 
relaxed  using  the  QM/MM  method  with  30  meV/A  and  1  meV  as  the  force  and  energy  convergence 
criterion  respectively.  The  QM  region  had  196  atoms  containing  the  dislocation  core  with/without 
an  oxygen  impurity  and  the  MM  region  consisted  of  1 1 ,926  atoms  capturing  the  long-range  elastic 
field  of  the  dislocation.  The  overall  dimensions  of  different  regions  were:  30  A  x  27  A  x  9.26  A 
for  the  QM  region  and  143.67  A  x  155.07  A  x  9.26  A  for  the  MM  region.  The  simulations  were 
performed  on  a  rectangular  cell  having  x,  y,  and  z-axes  oriented  along  the 
[1210],  [1010],  and  [0001]  directions,  respectively.  The  length  along  the  dislocation  line  [0001] 
was  set  to  2c  ( c  =  4.645  A)  to  avoid  the  interaction  of  oxygen  atoms  due  to  periodic  boundary 
conditions. 

For  the  DFT  part,  projector  augmented  wave  (PAW)  [21]  potentials  were  used  to  represent  the 
nuclei  core  with  valence  electrons  on  5  and  d  orbitals  for  a-Ti  and  valence  electrons  on  5  and  p 
orbitals  for  oxygen  atoms.  Exchange  and  correlation  was  treated  with  GGA  using  the  PBE  [22] 
form  with  an  energy  cutoff  of  289  eV  and  the  Monkhorst  Pack  k-point  mesh  of  1  x  1  x  5  along  the 
[1210]  ,  [1010],  and  [0001],  respectively.  For  GSFE  calculations  the  k-point  mesh  was  18  x  1  x 
10  along  the  [1210],  [1010],  and  [0001],  respectively.  A  single  oxygen  atom  was  introduced  on 
the  prismatic  shearing  plane  with  1/6  monolayer  (ML)  oxygen  concentration;  the  slip  direction  is 
along  [1210].  The  oxygen  atom  was  free  to  relax  in  all  directions;  whereas,  Ti  atoms  were  free  to 
move  along  the  [1010]  direction  only. 
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Figure  1 :  C-NEB  path  and  energy  barriers 
for  an  oxygen  atom  diffusing  in  bulk  a-Ti 
between  the  interstitial  sites. 


Reaction  Co-ordinates 


octahedral 


First  we  focus  on  oxygen  diffusion  in  bulk  a-Ti.  There  are  three  interstitial  positions, 
octahedral,  crowdion,  and  hexahedral  sites,  which  the  O  may  occupy,  and  the  preferred  position  is 
the  octahedral  site  [8],  The  hexahedral  site  is  on  the  basal  plane;  whereas,  the  octahedral  and 
crowdion  sites  are  between  the  A  and  B  stacking  along  the  [0001]  direction  [8],  The  transition 
pathways  and  energy  barriers  between  the  different  interstitial  sites  for  a  diffusing  O  atom  have 
been  determined  using  the  C-NEB  method  with  4  intermediate  images.  Figure  1  shows  the  energy 
barriers  for  the  O  atom  diffusing  from  an  octahedral  site  to  an  adjacent  crowdion  site  (O-C)  and 
subsequently  from  a  crowdion  position  to  a  hexahedral  position  (C-H)  and  finally  from  a 
hexahedral  position  to  an  octahedral  position  (H-O).  Also  in  Figure  1,  the  diffusion  pathway  is 
illustrated  by  showing  the  oxygen  positions  along  the  transition  path.  The  potential  energy  of  the 
hexahedral  site  is  1 .43  eV  higher  than  that  of  the  octahedral  site.  Moreover  the  excess  energy  of 
the  crowdion  site  over  the  octahedral  site  is  1.95  eV.  Both  crowdion  and  hexahedral  sites  are  local 
energy  minima  and  thus  excess  thermal  energy  is  required  for  the  O  atom  to  reach  the  octahedral 
sites.  The  diffusion  energy  barriers  from  the  crowdion  and  the  hexahedral  sites  to  the  octahedral 
sites  are  120  meV  and  -550  meV  respectively.  Both  compare  well  with  the  work  of  Wu  and  Trinkle 


[8], 
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Figure  2:  a)  Oxygen  binding  at  the 
prismatic  dislocation  core  with  the  plane  of 
sight  along  the  [0001]  direction,  b) 
Oxygen  binds  with  six  Ti  atoms  at  the 
dislocation  core.  Atoms  are  marked  by 
number  on  the  basal  plane. 


Next  we  study  oxygen-dislocation  interaction  in  a-Ti.  To  this  end  we  first  need  to  determine 
the  relative  site  preference  of  oxygen  at  the  dislocation  core.  The  total  energy  of  a  single  oxygen 
atom  placed  at  various  interstitial  sites  near  and  away  from  the  dislocation  core  has  been  calculated. 
Figure  2a  shows  the  optimized  equilibrium  geometry  of  the  oxygen  at  the  prismatic  edge 
dislocation  in  a-Ti  viewed  along  the  [0001]  direction.  It  is  found  that  the  local  atomic 
configuration  changes  from  an  HCP  octahedral  to  a  BCC  octahedral  structure.  The  oxygen  atom 
binds  with  four  Ti  atoms  on  the  same  basal  plane  and  with  two  Ti  atoms  above  and  below  the 
plane.  The  O  position  is  the  same  as  in  the  hexahedral  site  in  bulk  a-Ti  except  that  now  there  are 
four  atoms  on  the  basal  plane  (see  Figure  2b)  instead  of  three  (Figure  1). 


Figure  3:  a)  Diffusion  path  of  an  oxygen  atom  along  the  dislocation  line  ([0001]  direction)  of  the 
prismatic  edge  dislocation  with  ABAB  stacking  sequence,  b)  The  C-NEB  path  for  an  oxygen  atom 
diffusing  along  the  dislocation  line. 
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Having  obtained  the  stable  oxygen  site  at  the  dislocation  core,  we  proceed  to  determine  the 
oxygen  diffusion  pathway  and  energy  barriers  along  the  dislocation  core  using  the  C-NEB  method 
with  16  intermediate  images.  Figure  3a  shows  the  initial  path  for  an  oxygen  atom  to  diffuse  over 
a  distance  of  c  (4.645  A)  in  the  dislocation  line  direction  of  [0001]  with  ABAB  stacking.  The 
energy  barrier  and  diffusion  path  for  the  pipe  diffusion  of  oxygen  along  the  dislocation  line  are 
depicted  in  Figure  3b.  Also  the  migration  pathway  is  illustrated  by  showing  the  oxygen  position 
along  the  path  in  Figure  3b.  The  energy  barrier  for  the  diffusion  path  from  one  BCC  octahedral 
position  to  another  BCC  octahedral  position  is  1.24  eV  which  is  much  smaller  than  the  bulk 
diffusion  barrier  of  2.0  eV.  To  put  the  results  in  perspective  we  can  estimate  the  required 
temperature  for  an  oxygen  atom  jumping  from  one  octahedral  site  to  the  next  in  the  bulk  a-Ti  and 
for  the  pipe  diffusion  along  the  dislocation  line.  Assuming  the  phonon  frequency  as  1013  s"1  and 
for  a  unit  probability,  the  required  temperature  for  the  bulk  diffusion  is  775  K,  much  higher  than 
that  of  the  pipe  diffusion  at  496  K. 


Figure  4:  a)  GSFE  as  a  function  of  shear  displacement  along  the  [1210]  direction  of  a-Ti  with 
and  without  1/6  ME  oxygen  on  the  slip  plane,  b)  A  3D  iso-surface  plot  for  charge  density  showing 
the  effect  of  O  at  the  prismatic  stacking  fault.  The  orange  and  black  atoms  represent  oxygen  and 
Ti  atoms  respectively.  The  yellow  and  cyan  iso-surfaces  represent  charge  accumulation  and 
depletion  respectively. 
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Next  we  examine  the  effect  of  oxygen  on  the  deformation  behavior  a-Ti  using  the  SVPN 
model  in  conjunction  with  the  ab  initio  determined  GSFE.  Figure  4a  illustrates  the  GSFE  curves 
as  a  function  of  the  shear  displacement  in  the  [1210]  direction  with  (the  oxygen  concentration  is 
1/6  ML)  and  without  the  interstitial  oxygen.  The  unstable  stacking  fault  energy  is  increased  from 
237  mJ/m2  to  303  mJ/m2  (28%  increase)  in  the  presence  of  oxygen  on  the  shear  plane.  The 
dislocation  core  width  and  the  Peierls  stress  were  calculated  using  the  SVPN  method  with  ab  initio 
GSFE.  The  Peierls  stress  for  the  prismatic  edge  dislocation  with  and  without  oxygen  was  found  to 
be  97  MPa  and  24  MPa  respectively  (-400%  increase).  The  simulations  also  revealed  that  the 
prismatic  dislocation  becomes  narrower  in  the  presence  of  oxygen,  from  1.1  A  to  0.9  A.  This 
narrowing  of  the  core  width  is  consistent  with  the  increase  of  the  Peierls  stress  since  the  Peierls 
stress  depends  exponentially  on  the  ratio  between  the  dislocation  core  width  and  the  atomic 
spacing  along  the  dislocation  line. 

The  mechanical  properties  of  a  material  are  determined  by  the  nature  of  atomic  bonding; 
therein,  both  the  shear  strength  and  cohesive  energy  are  affected  by  the  directionality  and  strength 
of  chemical  bonds.  Therefore,  understanding  the  change  in  charge  density  resulting  from  a  solute 
atom  can  shed  light  on  the  strengthening  mechanism  of  the  impurity  (see  [23]).  Figure  4b  presents 
a  3D  iso-surface  plot  for  the  differential  charge  density  in  the  presence  of  the  oxygen  atom.  The 
yellow  and  cyan  iso-surfaces  represent  charge  accumulation  and  depletion  respectively.  Compared 
to  the  charge  density  distribution  in  pure  a-Ti  the  main  change  due  to  oxygen  happens  around  the 
octahedral  site  as  seen  in  Figure  4b.  The  charge  densities  between  Ti  atoms  were  evidently 
increased  with  the  addition  of  oxygen  resulting  in  an  increase  in  the  shear  modulus  which  affects 
the  dislocation  mobility. 

In  summary,  the  oxygen-metal  interaction  in  a-Ti  was  examined  including  oxygen  diffusion 
barriers  and  oxygen  effects  on  the  prismatic  edge  dislocation  ({1010}(1210))  motion.  It  was  found 
that  oxygen  prefers  to  occupy  the  BCC  octahedral  site  at  the  dislocation  core.  Moreover,  the 
oxygen  diffusion  barrier  along  the  edge  dislocation  line  from  a  BCC  octahedral  position  to  an  HCP 
octahedral  position  is  1 .24  eV  which  is  much  smaller  than  the  bulk  diffusion  barrier  of  2.0  eV.  We 
show  that  a  one-sixth  monolayer  of  oxygen  addition  increases  the  Peierls  stress  by  four-times  and 
reduces  the  dislocation  core  width  by  18%.  The  charge  densities  between  Ti  atoms  were  increased 
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with  the  addition  of  oxygen  at  the  prismatic  fault  resulting  in  an  increase  in  the  shear  modulus. 
The  calculated  hardening  effect  of  oxygen  and  oxygen  diffusion  barriers  are  consistent  with  those 
experimentally  observed  [3],  Overall  our  study  provides  critical  knowledge  towards  a 
comprehensive  understanding  of  oxygen  effects  on  the  deformation  behaviors  of  a-Ti. 

3.2  Grain  boundary  energy  and  structure  in  a-Ti 

Grain  boundaries  are  two-dimensional  defects  which  are  characterized  by  five  macroscopic 
degrees  of  freedom  [24-31].  Saylor  et  al.  [24]  analyzed  the  GB  character  distribution  (GBCD)  in 
commercially  pure  Al  and  suggested  that  boundaries  with  low  index  planes  occur  with  particularly 
high  frequency  in  polycrystals.  The  GB  structure-energy  correlation  in  several  FCC  metals  was 
investigated  by  Wolf  [32]  who  established  a  linear  relationship  between  GB  energy  and  volume 
expansion  per  unit  area  (grain  boundary  free  volume).  The  role  of  the  GB  plane  in  determining  GB 
energy  was  investigated  by  Holm  et  al.  [26],  who  performed  extensive  calculations  of  GB  energies 
in  Ni  and  Al.  Rohrer  [33]  established  that  the  coincidence  lattice  site  density  (I)  plays  only  a  minor 
role  as  a  determinant  of  GB  energy  in  several  systems  and  that  the  energies  of  the  crystal  surfaces 
that  meet  at  a  GB  is  a  better  indication  of  GB  energy.  Solanki  et  al.  found  that  certain  defects  are 
favored  at  a-Fe  grain  boundaries  and  that  these  species  affect  the  grain  boundary  cohesive  strength 
[34],  Experimentally,  grain  boundary  structure  has  been  observed  using  field  ion  microscopy  and 
high  resolution  transmission  electron  microscopy  [35-40].  The  grain  boundary  energies  can  be 
computed  through  theoretical  formulations  and  computational  methods.  Thus,  understanding  the 
atomic  structure  at  the  grain  boundary  can  provide  insight  into  the  grain  boundary  strength  as  well 
as  various  grain  boundary  dependent  phenomena,  such  as  dislocation  pileup  and  crack  initiation 
site  [41], 

Understanding  the  structure  and  energy  of  the  grain  boundary  system  is  crucial  for  engineering 
materials  intended  for  advanced  applications  because  grain  boundary  properties  can  vary  widely 
(coherent  twin  versus  low-angle  versus  high- angle  grain  boundaries).  Grain  boundaries  with 
certain  misorientation  angles  (and  typically  a  low  E  value)  correspond  to  “favored”  structural  units, 
while  other  boundaries  are  characterized  by  structural  units  from  the  two  neighboring  favored 
boundaries.  Also,  different  grain  boundaries  can  be  characterized  according  to  their  energy  and 
can  be  correlated  to  the  behavior  at  the  atomic  level.  In  the  first  year  of  this  project,  a  range  of 
grain  boundary  structures  and  energies  that  are  representative  of  some  of  the  variation  observed  in 
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the  grain  boundary  character  distribution  of  poly  crystalline  Ti  was  generated  to  be  used  in 
investigating  the  role  of  grain  boundary  character  on 


1 .  slip  transfer  behavior  (crack  initiation) 

2.  cohesive  zone  properties  (crack  growth) 

3.  non-Schmid  behavior  (dislocation  versus  twin,  crystal  plasticity  formulation) 

Initially,  a  database  of  61  minimum  energy  symmetric  tilt  grain  boundaries  (STGBs)  of  Ti, 
with  the  tilt  axes  as  [1210]  and  [0110]  was  generated  with  molecular  static  simulations,  which 
were  performed  using  the  classical  molecular  dynamics  code  Large-scale  Atomic/Molecular 
Massively  Parallel  Simulator  (LAMMPS)  [42],  Here,  the  analysis  cell  consisted  of  a  standard 
bicrystal  cell  with  a  single  grain  boundary  that  divides  the  HCP  crystal  into  two  single  crystals,  as 
shown  in  Figure  la.  The  initial  single  crystals  were  created  with  x,  y  and  z  along  the  [0110], 
[0001],  and  [1210]  directions,  respectively,  for  the  [1210]  tilt  axis;  and  along  the  [1210],  [0001], 
and  [0110]  directions,  respectively,  for  the  [0110]  tilt  axis.  Then,  the  upper  half  crystal  was 
rotated  clockwise  and  the  lower  half  crystal  counter-clockwise  by  angle  0  with  respect  to  the  tilt 
axis,  as  shown  in  Figure  5b.  Several  successive  rigid  body  translations,  followed  by  an  atom- 
deletion  technique  and  energy  minimization  using  a  non-linear  conjugate  method  [43 — 45]  were 
used  to  generate  the  final  relaxed  minimum  structure  with  the  grain  boundary  plane  along  the  x-z 
plane,  as  shown  in  Figure  5c.  This  procedure  was  replicated  to  generate  several  0  K  minimum- 
energy  grain  boundary  structures  for  Ti.  Here,  we  used  the  MEAM  potential  of  Hennig  et  al.  [46], 


v[0001] 


t  a)  b)  c) 


z[1210] 


x[0110] 


Figure  5:  (a)  Single  crystal  model  with  x,  y  and  z  along  the  [0110],  [0001],  and  [1210]  directions, 
respectively;  (b)  Upper  half  of  crystal  rotated  0  (-42°)  clockwise  and  lower  half  of  crystal  rotated 
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counter-clockwise;  and  (c)  Final  grain  boundary  structure:  rotated  crystal  is  minimized  with  1  pN 
force  on  each  atom  using  a  conjugate  gradient  algorithm.  Note:  TB  notes  the  twin  boundary. 

Figure  6  shows  both  the  grain  boundary  energy  as  a  function  of  misorientation  angle  for  the 
[1210]  STGB  system  and  the  variation  of  GB  energy  as  a  function  of  GB  geometry  in  a 
stereographic  triangle  representation,  which  is  widely  used  to  represent  cubic  systems  [32],  The 
trend  observed  for  the  grain  boundary  energy  as  a  function  of  a  misorientation  angle  is  comparable 
to  what  has  been  previously  reported  for  Ti  [47].  The  stereographic  representation  of  GB  energies 
for  the  entire  misorientation  range  helps  to  identify  GBs  that  show  a  local  minima  in  GB  energies 
in  the  energy-misorientation  plot  (refer  to  Figure  6b). 
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Figure  6:  (a)  Grain  boundary  energy  for  the  [1210]  symmetric  tilt  grain  boundary  system,  (b) 
Contour  plot  of  GB  energies  for  the  three  symmetric  tilt  systems  of  Ti  represented  using  polar  and 
azimuthal  angles.  The  polar  and  azimuthal  angles  correspond  to  the  degrees  of  freedom. 


The  structure-energy  correlation  can  provide  more  details  about  the  variation  in  grain  boundary 
energies,  as  each  grain  boundary  has  characteristic  structural  units  (SUs)  describing  its  atomistic 
morphology.  Low-angle  boundaries  can  be  represented  by  an  array  of  discrete  dislocations. 
However,  at  higher  misorientation  angles  (high-angle  grain  boundaries),  the  dislocation  cores 
overlap,  and  dislocations  rearrange  to  minimize  the  boundary  energy.  The  resulting  grain  boundary 
structures  are  often  characterized  by  grain  boundary  dislocations  or  SUs  [47 — 49].  Figure  7  shows 
the  spatial  distribution  of  atomic  excess  volume  for  (1011)[1210],  (2021)[1210], 
(2112)[0110],  and  (2116)[0110]  grain  boundaries.  Notice  that  the  atoms  far  away  from  the 
boundary  are  green,  indicating  that  there  is  no  atomic  volume  difference  over  the  bulk  lattice.  The 
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bulk  Voronoi  volume  was  found  to  be  17.57  A3  for  the  Ti.  The  excess  Voronoi  volume  is 
highest/lowest  (tensile  versus  compressive)  at  the  grain  boundary  center  and  converges  to  the  bulk 
Voronoi  volume  as  distance  from  the  grain  boundary  increases.  This  grain-boundary  metric  can  in 
turn  be  correlated  to  other  energetics  associated  with  the  grain  boundaries  to  derive  a  structure- 
property  relationship.  The  developed  grain  boundary  database  will  be  applied  to  investigating  the 
role  of  grain  boundary  character  on  (a)  slip  transfer  behavior  (crack  initiation);  (b)  cohesive  zone 
properties  (crack  growth);  and  (c)  non-Schmid  behavior  (dislocation  versus  twin,  crystal  plasticity 
formulation) 


p^  m  — 

a)  *4«  (1011)  *4***< 

/j  **  **  **  **  **  **  ** 

4*4  *4*4 *4*4*4* 4 * 
»*•*•*•*•*•*•*•*< 
*  **  **  **  M  *4  **  *•  * 

>4*4*4*4*4*4*4*44 
**************** 
>**  **  **********  *4 
4  **  *4  **  *4  *4  *4  **  *< 

l*  ************  **J 

4  *************** 

1.513(A3)|| 


-1.028(A3) 


b)  Sjl  (2021)  £?#/•« 

?4  4/fj*  ?*  *  *  4)f*  4  ^*  *  4 

,»4Tj**TT*4fJ*«T*4i 

**J***J**4*Jl**J**J 

jv*  35*4 

44  *  **  *  *  **  4  j  j*  4V,*  * 
A*  *££*  *^£*  *<lg&j4  *H*J 


r^*  4  J  J*  JjV  4  * V  *■ 


Figure  7:  Atomic  representation  of  (a)  (1011)[1210],  (b)  (2021)[1210],  (c)  (2 1 12)  [0110]  and 
(d)  (2116)[0110]  grain  boundaries  depicting  the  distribution  of  the  excess  Voronoi  volume.  The 
bulk  Voronoi  volume  was  found  to  be  17.57.  The  black  line  indicates  the  twin  boundary  plane. 


3.3  Atomistic  study  of  the  dislocation-grain  boundary  interaction 

The  short  fatigue  crack  propagation  behavior  is  significantly  altered  in  the  presence  of 
microstructural  features.  Grain  boundaries  play  a  fundamental  role  in  determining  the 
strengthening  mechanisms  as  they  act  as  an  impediment  to  dislocation  motion  [50,51],  The 
presence  of  grain  boundary  that  presents  a  large  barrier  to  plastic  flow  leads  to  generation  of 
dislocation  pileups.  A  simple  geometry  of  a  grain  with  the  basal  plane  inclined  to  the  applied  axis 
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and  a  TB-dislocation  interaction  model  are  shown  in  Figures  8a  and  8b,  respectively.  There  are 
four  possible  outcomes  of  slip-grain  boundary  interactions  (Figure  9)  i.e.  a)  a  direct  transmission; 
b)  a  direct  transmission  with  a  residual  dislocation  along  the  grain  boundary;  c)  an  indirect 
transmission  with  a  residual  dislocation  occurring  because  the  incoming  and  outgoing  slip  planes 
do  not  intersect;  and  d)  no  transmission  because  the  dislocation  is  absorbed  at  the  grain  boundary. 
Furthermore,  the  slip-grain  boundary  interactions  can  be  summarized  in  terms  of  the  Burgers 
vector  of  the  incident  ( bi ),  transmitted  (bt)  and  the  residual  dislocations  ( br ). 

*  bf  T  br 


Figure  8:  (a)  A  schematic  illustrating  dislocation  pile-up,  interaction  and  extraction  at  the  GBs 
that  leads  to  static  extrusions  across  the  boundary  in  the  form  of  ledges  and  steps,  and  (b)  a 
schematic  illustration  of  a  twin  lamella  with  the  dislocation  interaction  that  leads  to  static 
extrusions  across  the  boundary  in  the  form  of  ledges  and  steps. 


Figure  9:  Slip-grain  boundary  interaction  cases:  (a)  a  direct  transmission  and  dislocation  cross¬ 
slip;  (b)  a  direct  transmission  with  a  residual  dislocation  along  the  grain  boundary;  (c)  an  indirect 
transmission  with  a  residual  dislocation,  where  slip  planes  in  each  grain  do  not  intersect;  and  (d) 
no  transmission. 
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Several  criteria  [52-56]  have  been  developed  to  predict  the  outcome  for  the  DGB  interaction.  Shen 
et  al.  [54-56]  devised  a  set  of  rules  to  predict  the  slip  transmission  when:  a)  the  angle  of 
intersection  between  both  the  incoming  and  outgoing  slip  plane  is  minimized  and  b)  the  resolved 
shear  stress  acting  along  the  slip  direction  is  maximized.  Further,  based  on  in-situ  TEM 
experiments,  Lee  et  al.  [53]  proposed  an  additional  criterion  to  the  model  of  Shen  et  al.  [54-56] 
that  ensures  the  outcome  of  slip  transmission  is  based  on  the  minimum  residual  dislocation  along 
the  GB,  also  known  as  the  Lee-Robertson-Bimbaum  (LRB)  criteria.  The  geometric  condition  for 
the  slip  transmission  can  be  expressed  by  the  relative  orientations  of  the  slip  planes  and  the  GB 
plane  normal  in  the  following  manner. 

M  —  lin.  lout  (1) 

where,  l  is  the  unit  normal  at  the  intersection  of  the  incoming  and  outgoing  slip  planes  with  the 
GB  plane. 


Figure  10:  a)  slip  band  cracks  across  the  grain  boundary  in  the  case  of  grain  16  and  grain  17  and 
b)  deformation  twins  in  grain  2  appear  to  have  been  nucleated  at  the  grain  boundary  with  grain  1 . 
The  slip  bands  in  grain  1  and  the  deformation  twins  in  grain  2  are  correlated.  This  was  adopted 
from  Wang  et  al.  [57]. 

Experimental  verifications  of  these  models  are  difficult  because  of  microstructural  heterogeneities 
introduced  during  processing  which,  along  with  dislocation  pile-ups,  can  significantly  alter  the 
DGB  interaction.  Wang  et  al.  [57]  showed  formation  of  cracks  due  to  the  interaction  of  twin-twin 
and  slip  band  formation  as  seen  in  Figure  10a.  Moreover,  he  also  has  shown  deformation  twins 
nucleation  from  the  grain  1 -grain  2  grain  boundary  (Figure  10b)  and  propagates  towards  the  grain 
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boundary  of  grain  2-grain  3  (Figure  10b),  based  on  the  difference  in  the  thickness  of  the  twin  at 
the  grain  boundary  of  grain  1 -grain  2  as  compared  with  the  grain  boundary  of  grain  2-grain  3. 
Furthermore,  a  slip-twin  correlation  with  the  help  of  the  Luster-Morris  parameter  [58];  wherein, 
the  prismatic  slip  interacts  with  the  grain  boundary  and  nucleates  a  deformation  twin  in  the 
adjacent  grain  even  though  the  c-axis  of  the  grain  is  oriented  at  45°  to  the  loading  axis  which  is 
unfavorable  for  twinning  according  to  the  Schmidt  factor.  Our  main  objective  was  to  quantify  all 
of  the  possible  outcomes  of  slip/  twin  interaction  with  the  grain  boundary  interaction. 


Further,  the  atomic  level  details  related  to  the  interaction  of  individual  grain  boundaries  with 
dislocations  still  remain  unclear  due  to  insufficient  resolution.  For  example,  the  GB  atomic 
structures  vary  spatially  which  could  change  the  pining  strength  of  impeding  dislocations.  This  has 
motivated  several  atomistic  studies  in  BCC  [59,60],  FCC  [61,62],  and  HCP  [63,64]  metals  to 
validate  these  analytical  models.  In  particular,  previous  atomistic  [61,65]  studies  have  shown  that 
slip  transmission  was  observed  to  violate  the  LRB  criteria  by  transmitting  on  slip  planes  that  did 
not  have  the  highest  resolved  shear  stress.  Based  on  these  findings,  it  was  suggested  that  the 
underlying  atomic  structure  at  the  DGB  interaction  site  plays  a  key  role  in  determining  the  outcome 
[61,62,65].  Despite  the  large  body  of  work  on  DGB  interactions,  there  is  still  a  lack  of  systematic 
study  on  the  role  of  grain  boundary  structure  especially  in  a-Ti. 

Here,  we  carried  out  atomistic  simulations  to  examine  the  DGB  interactions  across  several  grain 
boundaries  in  a-Ti  using  interatomic  potential.  The  energy  barrier  for  slip  transmission  across  the 
GB  was  used  to  quantified  by  defining  a  control  volume  at  the  site  of  the  DGB  interaction  [62], 
The  defected  atoms  (centrosymmetry  parameter  [66]  >  0.5)  within  the  control  volume  were  used 
to  estimate  the  energy  barrier  for  slip  transmission.  The  evolution  of  the  net  change  in  energy  of 
the  defected  atoms  during  the  loading  process  was  quantified  by  comparing  instantaneous  energy 
(Ej-)  with  reference/initial  atomic  energy  (Er).  The  net  energy  was  normalized  with  the  atomic 
volume  occupied  by  defect  atoms  in  the  reference  configuration: 


(2) 


where  V  represents  the  total  volume  of  the  defected  atoms  within  the  control  volume  in  the 
reference  configuration. 
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Figure  lira)  The  dislocation  transmission  energy  barrier  across  a  grain  boundary  as  a  function  of 
the  grain  boundary  energy  in  Titanium,  b)  The  atoms  were  colored  according  to  the  local  atomic 
shear  strain  (yXy)  on  a  scale  of  0  to  0.4  showing  trace  of  dislocation  glide  and  nucleation  of  a  twin. 


In  this  work,  seven  different  grain  boundaries  from  the  [2110]  and  [0110]  tilt  axes  were  deformed 
to  study  the  DGB  interaction.  A  strong  inverse  relationship  between  the  energy  barrier  for  slip 
transmission  ( Ebarrier )  and  the  initial  GB  energy  (EGB)  was  found  Figure  1 1.  In  other  words,  GBs 
with  lower  interfacial  energy  demonstrated  a  higher  barrier  for  slip  transmission,  which  is  in 
agreement  with  previous  study  [62],  Figure  1  lb  shows  the  interaction  of  the  Basal  dislocation  in 
grain  1  (soft  grain)  with  the  grain  boundary.  Subsequently,  twin  nucleates  from  the  DGB 
interaction  site  and  propagates  in  neighboring  grain  (grain  2). 


The  energy  barrier  calculations  for  the  dislocation  nucleation  from  various  GBs  and  subsequent 
slip  transfer  can  be  cast  into  a  crack  initiation  criterion  (Cl),  which  is  a  function  of  the  slip  band 
geometry,  Ebarr ier,  and  slip  increment.  The  total  energy  for  the  slip  configuration  is  given  by 

E  xj)bibjh(p  —  Pint)  ^barrier,  (1) 

where  (xj  —  x;)  are  the  dislocation  positions,  b  is  the  Burgers  vector,  bibjh.  is  the  slip  band 
geometry,  a  is  the  model  dimension  constant,  and  (p  —  pint)  is  the  dislocation  density  term.  The 
crack  initiation  is  predicted  by  taking  the  first  derivative  of  the  total  energy  with  respect  to  the 
dislocation  increment  as 

Cl  =  f  (2) 

dx 
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To  further  understand  the  slip  transfer  and  deformation  behavior  of  poly  crystalline  materials, 
we  need  to  investigate  the  mechanical  behavior  of  the  bulk,  i.e.,  the  grain  interiors  and  the  defects, 
which  in  our  case,  are  the  grain  boundaries.  Nanoindentation  offers  a  great  advantage  in  terms  of 
isolating  and  probing  a  small  volume  to  investigate  the  deformation  and  slip  transfer  across  the 
grain  boundaries.  The  load-displacement  data  from  nanoindentation  coupled  with  the  orientation 
information  obtained  through  EBSD  can  help  in  characterizing  a  large  set  of  grain  boundaries 
found  experimentally.  We  pick  commercially  pure  Titanium  (CP-Ti)  as  our  primary  material 
system.  Then  we  propose  to  introduce  solutes  (alpha  stabilizing  agents  such  as  O)  to  CP-Ti  through 
processing  techniques  since  it  is  known  that  the  deformation  is  very  sensitive  to  the  changes  in 
composition  of  the  alloy. 

This  experimental  study  is  primarily  composed  of  three  levels  per  composition  of  the  alloy. 
Firstly,  the  critical  resolved  shear  stress  values  for  slip  on  the  three  main  systems,  namely,  basal, 
prismatic,  and  pyramidal  are  analyzed  to  understand  the  plastic  deformation  of  CP-Ti.  The  critical 
resolved  shear  stress  can  be  estimated  using  the  load  at  the  pop-in  which  is  observed  during 
indentation  of  the  grain  interior  region.  The  Hertzian  theory  dictates  the  elastic  response  and  the 
elastic-plastic  transition  identified  with  a  discontinuity  in  the  load-displacement  curve.  This  is 
known  as  a  “pop-in”  event  [67].  The  pop-ins  are  characteristic  of  initiation  of  plastic  flow  or,  in 
other  words,  activation  of  dislocation  sources  in  the  indented  volume.  Secondly,  by  probing  the 
regions  very  close  to  the  grain  boundary,  an  additional  pop-in  event  can  be  observed  that  is 
indicative  of  the  interactions  between  the  grain  boundaries  and  the  dislocations  source.  This 
secondary  pop-in  occurs  at  higher  loads  and  thereby  after  appreciable  plastic  deformation  and 
corresponds  to  the  transmission  of  the  piled  up  dislocations  (after  the  initial  pop-in  event)  across 
the  grain  boundaries[68].  Previous  studies  have  identified  the  emission  of  dislocations  on  the  other 
side  of  the  boundary  for  low  angle  boundaries.  However,  for  high  angle  boundaries,  the 
transmission  was  not  observed  due  to  the  predominance  of  GB  sliding[69,70].  It  has  been  observed 
that  pop-in  occurs  for  a  critical  ratio  between  the  plastic  zone  size  beneath  the  indenter  to  the 
distance  between  the  indent  and  the  grain  boundary  plane  for  a  given  misorientation  angle  [71]. 
Furthermore,  the  stress  intensity  factor  for  slip  transfer  was  reported  to  increase  with  the  increase 
in  the  misalignment  of  slip  systems  (m’  parameter)[72].  It  has  been  pointed  out  by  Britton  et  al. 
[72]  that  it  is  necessary  to  identify  the  detailed  stress  state  of  the  indenter  and  correlate  the 

Page | 16 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


parameters  to  the  grain  boundary  character  better  to  completely  understand  the  process  of  slip 
transfer  during  nanoindentation.  Lastly,  pileup  topographies  will  also  be  imaged  through  atomic 
force  microscopy  to  identify  the  anisotropy  in  plastic  deformation  for  the  various  slip  systems.  By 
isolating  the  bulk  contribution  from  the  grain  boundary  contribution,  we  can  quantitatively 
determine  the  role  of  a  grain  boundary  on  the  slip  transfer  process. 


Displacement  zone 


Figure  12:  Schematic  illustration  of  (a)  load-displacement  curve  with  initial  and  grain  boundary 
pop-in,  and  (b)  direct  transmission  of  dislocations  occurring  across  a  grain  boundary  (adapted 
from  [72]) 

CP-Ti  samples  for  characterization  and  testing  were  machined  through  electrical  discharge 
machining,  and  then  mechanically  ground,  polished,  and  finished  with  electrochemical  polishing. 
As  far  as  the  composition  is  concerned,  the  oxygen  content  was  quoted  as  0.2  wt.%  by  the 
manufacturer.  The  samples  were  carefully  stored  to  minimize  the  effect  of  ambient  air  on  the 
surface.  As  the  next  step,  we  characterized  the  local  microstructure  of  the  sample  through  high 
resolution  EBSD  analysis  for  an  area  spanning  600  pm  x  400  pm  (Figure  13).  The  data  was  cleaned 
to  reduce  erroneous  points.  Subsequently,  regions  of  interest  were  identified  for  nanoindentation. 
Hysitron  Triboscope  were  used  to  perform  indents  on  the  sample  using  sphero-conical  and 
berkovich  tips.  To  ensure  the  accuracy  of  the  test,  a  sufficient  number  of  indents  were  made  in  the 
grain  interior.  Care  was  taken  to  reduce  the  effect  of  grain  boundaries  and  other  grains  for  step  1 
that  is,  computing  the  critical  resolved  shear  stress  of  the  grain  interior  (bulk).  The  other 
geometrical  factors  were  calculated  using  the  orientation  information  obtained  from  EBSD 
analysis.  An  illustration  indicating  the  grain  boundary  and  grain  indentation  scheme  is  represented 
in  Figure  13.  Topography  around  the  indents  were  measured  using  AFM’s  tapping  mode  on  a 
scanning  probe  microscope. 
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Figure  13:  (a)  High  resolution  EBSD  map  of  a  region  from  a  CP-Ti  sample,  and  (b)  representation 
of  grain  interior  (GI)  and  grain  boundary  (GB)  indentation  scheme 


Figure  14  shows  the  load  versus  displacement  curves  of  two  of  the  indents  made  at  the 
grain  interior  and  near  the  grain  boundary  for  CP-Ti.  The  indent  of  the  GB  was  made  2  pm  away 
from  the  grain  boundary.  Displacement  jumps/ strain  bursts  are  evident  from  the  curves  at  various 
loads.  The  first  jump  occurs  at  about  mN,  and  it  is  observed  for  all  the  indents  made  within  the 
grain  and  close  to  the  grain  boundary.  The  pop-in  is  characteristic  of  the  transition  from  the  elastic 
to  the  plastic  regime  which  can  be  associated  with  the  stress  required  for  dislocation  nucleation. 
Multiple  indentations  on  several  grains  were  used  to  estimate  the  critical  resolved  shear  stress  for 
various  slip  systems.  The  average  CRSS  values  for  three  slip  systems  (prismatic:  basal:  pyramidal 
<c+a>)  was  computed  to  be  181:209:474  MPa  which  is  approximately  1:1.26:2.45.  This  is 
consistent  with  the  values  reported  in  the  literature  (1:1.2:2.62)  where  prismatic  was  found  to  be 
easily  activated  for  slip  in  CP-Ti  material  system.  In  case  of  indents  close  to  the  grain  boundary, 
there  occurs  a  second  pop-in  similar  to  Figure  14  (b)  at  higher  loads  (in  this  case)  which  is 
dependent  on  the  distance  of  the  grain  boundary  to  the  indent  and  orientation  of  the  GB.  The  radius 
of  the  elasto-plastic  boundary  was  estimated  using  the  equation  proposed  by  Kramer  which  takes 
into  account  the  Pop-in  load,  and  the  yield  strength  of  the  material.  Several  values  of  c/d  were 
estimated  from  the  indentation  data  which  suggests  a  critical  value  of  0.2  to  yield  pop-ins  (based 
on  Johnson’s  cavity  model).  The  critical  stress  intensity  factor  for  slip  transfer  was  computed  based 
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on  the  formulation  of  Wang  et  al.  suggests  a  range  of  values  to  facilitate  slip  transfer  between 
grains  of  varying  orientation  (0.003  and  0.41  MPam0  5).  Although  a  range  of  values  are  reported, 
these  correspond  to  the  different  GBs  in  the  material  system.  The  slip  systems  which  are  highly 
aligned  facilitate  slip  transfer  and  exhibit  pop-in  behavior  have  a  M  >  0.98. 


to 

8 

b 

4 

2 

0  200  400  0 

Displacement  (nm) 

Figure  14:  Load-displacement  curves  highlighting  the  pop-ins  in  (a)  grain  interior  (GI)  and  (b) 
near  a  grain  boundary  (GB)  in  CP-Ti  Grade  2  material. 

In  case  of  Ti-0  sample,  several  levels  of  O  was  introduced  through  boost  diffusion 
procedure  where  a  CP-Ti  sample  was  subjected  to  heat  treatment  under  vacuum  conditions.  Prior 
to  indentation,  an  EBSD  analysis  was  made  to  identify  the  texture  and  nature  of  microstructure  of 
the  sample  (Figure  15).  Multiple  indentations  on  the  grain  interior  and  near  the  boundary  were 
made  and  the  first  pop-in  loads  for  indents  were  observed  to  increase  with  the  increasing  O  content 
(Figure  15  (b)).  Further,  the  hardening  effect  of  the  O  solute  is  evident  from  the  load-displacement 
curves  and  the  absence  of  twinning  indicates  the  increase  in  the  twin  formation  energy  due  to  the 
presence  of  O.  The  CRSS  values  for  the  slip  systems  revealed  a  similar  ratio  as  that  of  pure  Ti 
sample  where  pyramidal  was  found  to  be  less  activated  than  prismatic  and  basal  slip 
systems(pyramidal<bsasl<prismatic).  The  critical  stress-intensity  factor  for  slip  transfer  revealed 
a  range  of  values  between  0.001  and  0.1  MPam0  5  for  the  grain  boundaries  indicating  that  it  is 
easier  for  dislocation  motion  compared  to  nucleation  in  the  presence  of  O  consistent  with  atomistic 
simulations.  In-situ  SEM  characterization  of  one  of  the  indents  indicated  in  Figure  15  (c)  of  Ti-0 
sample  reveals  the  pile-up  of  the  material  around  the  indent. 
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Figure  15:  (a)  High  resolution  EBSD  map  of  a  region  from  a  Ti-0  sample,  (b)  variation  in  pop- 
in  loads  in  Ti-0  with  respect  to  O  concentration,  and  (c)  pile  up  of  material  around  the  indent  for 
in- situ  indentation  of  a  grain  interior. 

3.4  Extended  finite  element  framework  for  the  small  crack  growth 

The  material  fatigue  strength  of  engineering  alloys  is  significantly  affected  by  microstructural 
variation.  For  example,  fatigue  crack  initiation  and  early  growth  are  primary  contributing  factors 
to  the  alloy’s  total  lifespan,  and  these  mechanisms  strongly  depend  on  grain  boundaries  and  other 
defects.  Furthermore,  microstructural  variations  such  as  relative  grain  orientations  influence  the 
cohesive  behavior  of  the  grain  boundaries,  which  leads  to  the  crack  nucleation  and  propagation  in 
an  intergranular  and/or  transgranular  (see  Figure  16)  manner  [73-76], 


Figure  16:  Intergranular  and  transgranular  small  cracks  [73] 
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The  extended  finite  element  method  (XFEM)  provides  a  powerful  tool  for  finite  element 
modeling  of  interfaces  independently  from  the  computational  mesh  [77].  To  describe  the  role  of 
microstructure  in  material  failure,  some  methods  have  been  proposed  among  which  the  cohesive 
finite  element  method  (CFEM)  has  received  much  attention  [78],  Usually,  empirical  data  for 
surface  decohesion  is  employed  in  CFEM  to  effectively  simulate  crack  initiation  or  propagation  in 
material.  Nevertheless,  empirical  cohesive  laws  might  be  incapable  of  describing  decohesion 
behavior  accurately  at  sub-micron  scales  [78],  As  an  alternative,  atomistic  simulations  can  be  used 
to  characterize  the  cohesive  behavior  of  a  variety  of  grain  boundaries  in  order  to  be  used  in  finite 
element  simulations. 


3.4.1  XFEM  Formulation 

As  reported  before,  the  suitable  XFEM  displacement  field  for  modeling  a  microstructure  having 

strong  discontinuities  is: 
nn  ng 

u(x)  =  ^  ^  N/(x)Hi(x)ui/  (1) 

1=1  i= 1 

where  in  the  above  equation,  NG  is  the  number  of  grains  in  the  microstructure  and  Nn  is  is  the 

number  of  nodes.  Ht  is  a  Heaviside  function  representing  grain  i  and  is  defined  as  follows: 

x  >  0  (2) 

10  x  <  0 

For  some  values  of  I  and  i  and  for  all  values  of  x,  N,(x)//j  =  0.  For  these  combinations,  u</  is 
prescribed  to  be  zero  and  is  not  stored  in  the  global  degrees  of  freedom.  The  strong  form  of 
equilibrium  equations  and  prescribed  boundary  conditions  are: 


H 


■(x)  =  (J 


V.  a  =  0 

on  0 

(3) 

a  =  C:  e 

(4) 

a.  n  =  t f 

on  dQ^ 

(5) 

©  — 

3 

II 

s 

- 

on  dQu. 

(6) 

C  and  e  in  Eq.  4  are  the  rank  four  elasto-plastic  stiffness  tensor  and  strain  field,  respectively. 
Additionally,  at  the  grain  boundaries  of  each  grain  (dQgb),  we  require  that: 

<T.n  =  U  on  dQgbi  (7) 

for  i  =  1, 2,  ■■■  ,Ng.  We  assume  a  traction-separation  constitutive  law  at  the  grain  boundaries  of 
each  grain  to  have  the  following  form: 
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t i  =  Ti.  lutJ  (8) 

where  T  is  a  second-order  constitutive  tensor  and  l  ]  denotes  the  displacement  jump  at  the  grain 
boundary.  Displacement  jump  ([u])  in  Eq.  8  for  two  adjacent  grains  such  as  Grain  a  and  Grain  b 
at  point  a  (or  point  b )  on  the  interface  (grain  boundary  ab )  is  defined  as  follows  (see  Figure  17): 
ttuOOIafc  =  Ua(x)  -  Ub(x)  (3) 


Figure  17:  Discontinuous  grain  boundary 


The  weak  form  equation  is  simplified  as  follows: 


/ 


C:  €, 


dv  +  ^ 


I 


wa- T, 


ab 


■  (ua  -  u„)  ds 


-I 


Wa.t"  ds 


(4) 


Qa  ®  d^gba  dQa 

where  subscript  b  is  the  number  of  any  grain  sharing  an  interface  with  grain  a  and  eais  the  strain 
field  corresponding  to  the  virtual  displacement. 


3.4.2  Finite  element  discretization 

The  discretized  form  of  the  equilibrium  equation  can  be  derived  from  the  weak  form  in  Eq.  9  and 
the  enriched  displacement  field  of  Eq.  (1): 

fint  =  Mext  +  fcoh  (10) 

where: 


fint=  j  HaBT{C}BHadv.ua 

a  a 

fC°h  =  Z  /  HaNT ' Tab ' (N//flUfl  “  NHbUb)  ds 

b  dQgba 


(11) 

(12) 
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(13) 


In  the  above  equations,  fint,  fcoh,  and  fext  are  the  internal,  cohesive,  and  external  forces  in  the 
system,  respectively.  In  order  to  trace  the  equilibrium  path,  Eq.  9  should  have  zero  residual: 


R(ua,  A)  =  Afext  +  fcoh— fint  =  o 


(5) 


Eq.  14  has  to  be  solved  for  both  ua  and  A  simultaneously  using  an  iterative  non-linear  solver. 
Nevertheless,  as  previously  stated,  an  arc  length  method  is  required  to  simulate  the  response  of  the 
system  when  there  is  a  reduction  in  the  material  stiffness  since  displacement  control  or  load  control 
methods  alone  fail  to  converge  at  the  onset  of  global  softening.  Here  we  implemented  Crisfield 
and  Modified  Crisfield-Ramm  arc  length  methods  [79],  A  robust  arc  length  method  facilitates  the 
study  of  load-displacement  of  the  microstructure  until  a  full  fracture  occurs. 

3.4.3  Method  verification 

In  order  to  verify  the  method  implementation,  a  simple  three-grain  specimen  (Figure  18)  with 
normalized  dimensions  was  modeled  in  ABAQUS  and  compared  with  the  present  method.  In  this 
case,  a  consistent  linear  cohesive  law  at  the  grain  boundaries  was  used.  A  1000  ( N/m )  normal 
traction  was  applied  on  the  right  edge  of  the  specimen  while  the  left  edge  was  held  fixed.  High 
stiffness  without  softening  is  modeled  for  elements  in  contact  to  prevent  material  interpenetration 
at  the  grain  boundaries. 
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Figure  18:  Three-grain  specimen  with  normalized  dimension 
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Aluminum  6061  material  properties  were  chosen  for  the  bulk  material.  Contours  of  oxx  for  the 
two  solutions  are  depicted  in  Figure  19.  A  good  agreement  between  the  XFEM  and  ABAQUS 
results  is  observed.  To  achieve  a  smooth  strain  field  in  an  ABAQUS  solution,  the  nodes  of  the 
adjacent  elements  at  the  two  sides  of  any  interface  have  to  be  coinciding.  This  becomes  an 
extremely  burdensome  task  to  fulfill  in  order  to  mesh  more  complex  geometries.  Nevertheless,  this 
problem  is  eliminated  in  XFEM  by  implicitly  modeling  the  interfaces.  In  fact,  relatively  accurate 
results  could  be  achieved  with  an  even  fewer  number  of  elements  in  the  XFEM  solution. 
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Figure  19:  (a)  XFEM  solution  with  a  2500-element  structured  mesh,  and  (b)  ABAQUS 
solution  with  a  7500-element  unstructured  mesh 


3.4.4  Effect  of  local  grain  boundary  characteristic  on  the  crack  growth 

In  this  section  a  qualitative  parametric  study  is  performed  to  demonstrate  the  effect  of  different 
variations  of  grain  boundary  cohesive  behavior  on  small  crack  growth  (SCG)  in  polycrystalline 
material.  As  mentioned  before,  traction-separation  laws  follow  the  Needleman’s  model  [82]: 

16e2  A  \ 

/  (6) 


16e2  A 

T(A)  =  T0-^-—exp  I 


9  An 


where  T  is  the  traction,  A  is  the  separation  distance,  and  A0  is  the  maximum  separation  distance 
before  new  cracked  surfaces  are  formed  and  T0  is  the  maximum  traction  in  the  traction-separation 
curve.  Reference  energy  for  interface  decohesion  (separation  energy)  is  considered  to  be  G  =  0.85 
(J/m2)  with  Tq  =  4  GPa  and  A0=  4(A).  Here  a  microstructure  shown  in  Figure  20  is  simulated  for 
three  cases  of  different  combinations  of  grain  boundary  cohesive  zone  models  (Table  1). 
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Figure  20:  Microstructure  under  tensile  loading 

All  the  six  different  cohesive  zone  models  (CZM)  in  Table  1  are  considered  and  normalized  with 
respect  to  the  reference  cohesive  zone  model  (CZM  1)  and  depicted  in  Figure  21. 


Figure  21:  Different  normalized  cohesive  zone  models 


Table  1:  Cohesive  law  for  three  test  cases 


Grain  Boundary 

GB  12 

GB  14 

GB  25 

GB  24 

GB  34 

GB  35 

GB  45 
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Case  1 

CZM  2 

CZM  3 

CZM  1 

CZM  1 

CZM  5 

CZM  6 

CZM  1 

Case  2 

CZM  2 

CZM  2 

CZM  2 

CZM  1 

CZM  6 

CZM  1 

CZM  1 

Case  3 

CZM  4 

CZM  4 

CZM  4 

CZM  3 

CZM  1 

CZM  4 

CZM  4 

Grain  boundary  sliding  stiffness  is  taken  to  be  a  high  value  compared  to  the  normal  stiffness. 
Additionally,  high  penalty  stiffness  is  assumed  to  avoid  interpenetration  at  grain  boundaries. 
Because  the  intent  here  is  to  investigate  the  effect  of  different  cohesive  zone  models  on 
intergranular  fracture,  elastic  modulus  of  the  bulk  material  is  taken  to  be  very  high.  The  fully 
debonded  (crack)  surfaces  are  delineated  in  Figure  22  by  white  lines.  Uniform  load  is  being  applied 
on  the  right  edge  and  the  other  edges  are  constrained.  An  in  Figure  22  denotes  the  externally  applied 
load,  non-dimensionalized  by  T0,  and  s>,  is  the  average  strain  of  the  right  edge.  It  should  be  noted 
that  these  values  are  given  merely  for  comparison  reasons  and  they  do  not  represent  the  actual 
response  of  the  microstructure. 


Case  1 


(a)  a„=0.45  ,  e„=0.58e-4 


(b)  A„=0.52  ,  e„= 0.99e-4  (c)  2„=0.28  ,  e„=l .  12e-4  (d)  2„=0.07,  e„=2.06e-4 


Case  2 
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(a)  A„=0.24 ,  e„=0.27e-4  (b)  A„=0.14  ,  e„=2.46e-4  (c)  ).„={). 01 ,  e„=2.48e-4  (d)  /L„=0.05  ,  e„=2.49e-4 


Case  3 


Figure  22:  Effect  of  different  grain  boundary  CZMs  on  small  crack  nucleation  and  growth 

As  it  is  evident  from  Figure  22,  different  grain  boundary  properties  can  yield  different  patterns  for 
crack  nucleation  and  growth  in  microstructure.  This  indicates  the  importance  of  the  incorporation 
of  each  grain  boundary  characteristic  in  finite  element  simulations. 


3.4.5  Crystal  plasticity  formulation  including  twinning 

Incorporation  of  crystallographic  plastic  slips  in  the  crystalline  materials  has  been  long  established 
using  rigorous  phenomenological  constitutive  laws  to  capture  the  anisotropic  stress-strain 
response,  evolution  of  texture,  and  macroscopic  behavior  of  the  polycrystalline  materials  [1,  2], 
Originally  the  crystal  plasticity  constitutive  laws  were  developed  for  FCC  material.  Nevertheless, 
effect  of  deformation  twinning  mechanism  as  a  mode  of  plastic  deformation  in  many  materials 
such  as  Ti  cannot  be  simply  neglected.  In  order  to  take  into  account  the  effect  of  twinning,  a 
statistical  criterion  is  used  which  is  based  on  the  volume  fraction  of  twinned  zones  in  the  entire 
microstructure. 
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Current  configuration 


Figure  23:  Multiplicative  decomposition  of  the  deformation  gradient. 


The  crystal  plasticity  constitutive  law  used  in  this  study  is  based  on  multiplicative  decomposition 
(see  Figure  23)  of  total  deformation  gradient  where  F  =  FeFp,  where  Fe  and  Fv  are  the  elastic 
and  plastic  deformation  gradients.  The  time  rate  of  plastic  part  of  the  deformation  gradient  is 
calculated  by 

F  =  LPFP  (36) 


where 

L"  =  (i-  sf""  fh  zf‘p  rss-sm  +  zfwt"  fpytw,nso-twm  (”) 


is  the  velocity  gradient  of  the  plastic  deformation,  S£_slip  and  Sp_twin are  the  Schmid  tensors  for 
the  slip  system  a  and  twinning  system  /?,  respectively,  and  is  the  volume  fraction  of  and 
twinning  system  /?.  ytwin  is  a  constant  shear  strain  associated  with  twinning.  A  phenomenological 
flow  rule  is  used  to  calculate  the  rate  of  slip  on  slip  system  a  [3] 


Ya  =  Yo 


sign(ra ) 


(38) 


where  y0  is  a  reference  slip  rate,  m  is  the  rate  sensitivity  parameter,  and  ra  and  saare  resolved 
shear  stress  and  resistance  on  slip  system  a,  repectively.  Evolution  of  deformation  twin  volume 
fraction  is  obtained  by  a  power-law  relation 
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fp  = 


(39) 


where  and  are  resolved  shear  stress  and  twin  resistance  on  twin  system  /?,  repectively.  The 
evolution  of  the  slip  and  twin  resistance  in  the  above  formulations  are  updated  by 


NsliP  1/a'id  ,>fc 

<a  Y  )  Lk  Y 


(40) 


(41) 


Klip  and  Kwin  are  slip  and  twin  hardening  modulus,  respectively.  C,  a,  b  and  d  are  hardening 
parameters  that  control  the  contribution  of  different  resistance  mechanism  for  slip,  twin  and  twin- 
slip  interaction.  Further  details  on  the  physical  phenomenology  of  the  above  formulae  can  be  found 
in  the  literature  [3], 

3.4.6  Fatigue  test 

In  order  to  estimate  the  effect  of  oxygen  on  fatigue  life  of  Ti,  strain  controlled  fatigue  tests  were 
performed  on  125  mm  long  single  edge  notch  tensile  specimens  of  CP-Ti  and  Ti-0  with  an  R  ratio 
of  0.1  under  displacement  control  method  (max  strain  ~  0.1%).  The  notch  was  machined  using 
EDM  and  samples  were  polished.  Imaging  was  performed  using  a  microscope  camera  with  a  lOx 
magnification  to  identify  the  crack  initiation  and  growth  ahead  of  the  notch.  DIC  was  performed 
with  the  images  captured  post-deformation.  Microprobe  analysis  was  performed  post  deformation 
due  to  restrictions  of  sample  size. 
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Figure  24:  Single  edge  notch  tension  sample  geometry  used  in  fatigue  test  and  DIC  measurement 
scheme  where  images  were  captured  at  no  load  and  maximum  load  conditions. 

For  the  finite  element  simulations,  crystal  plasticity  parameters  were  first  calibrated  using  ab-initio 
calculations  on  dislocation  motion  energetics.  As  it  was  shown  in  the  previous  section,  presence 
of  oxygen  reduces  the  slip  activity  by  increasing  the  dislocation  motion  energy  on  some  of  the  slip 
systems  in  Ti,  including  the  dominant  prismatic  system.  Table  1  shows  how  oxygen  changes  the 
maximum  stacking  fault  energy  on  different  slip  systems  in  Ti  which  consequently  leads  to 
changes  in  the  corresponding  CRSS  values. 


Table  1:  Dislocation  motion  energy  (GSFE)  for  CP-Ti  and  Ti-O 


Slip  System 

Ti 

Ti-O  (20%  ml  oxygen  content) 

Ratio 

Prismatic 

237 

303.17 

1.28 

Basal 

278 

168 

0.60 

Pyramidal<a> 

761.21 

1026.7 

1.34 

Pyramidal<c+a> 

442 

430 

0.97 

This  can  significantly  contribute  to  the  increase  in  the  initial  yield  strength  of  Ti  in  presence  of 
oxygen,  which  has  been  depicted  in  Figure  25  and  compared  with  crystal  plasticity  FEM 
simulations. 
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Figure  25:  True  stress  vs.  true  strain  curves  for  CP-Ti  and  Ti-0  obtained  from  experiments  and 
compared  against  CPFEM  simulations. 

On  the  other  hand,  degrading  effect  of  oxygen  on  the  grain  boundaries  makes  them  preferable 
sites  for  nucleation  and  propagation  of  small  cracks  in  an  intergranular  fashion,  which  is 
consistent  with  the  experimental  results  (Figure  26). 


Figure  26:  DIC  images  and  underlying  microstructure  in  CP-Ti  (left)  and  Ti-0  (right).  The 
contours  indicate  magnitude  of  strain  during  fatigue  testing  along  the  loading  direction. 


CPFEM  simulations  were  carried  out  on  2D  notched  samples  with  the  calibrated  parameters. 
In  order  to  reduce  the  computational  costs,  the  process  zone  was  limited  to  a  region  close  to 
the  notch  and  the  rest  of  the  domain  was  considered  to  deform  in  an  isotropic  homogenized 
elasto-plastic  manner  (see  Figure  27). 
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Figure  27:  Computational  geometry  and  the  modeling  strategy  for  CPFEM. 

This  assumption  is  valid  since  the  notch  stress  concentration  effect  is  dominant  only  within  a  small 
region  ahead  of  notch  in  a  microstructure-sensitive  model.  This  robust  framework  can  capture  the 
complexity  of  a  real  microstructure,  effects  of  notch  size  and  characteristic  length  in  the 
distribution  of  slips  in  the  microstructure,  which  can  be  used  to  estimate  the  fatigue  life  of  the 
material.  The  normalized  total  plastic  shear  strain  contours  for  CP-Ti  and  Ti-0  tests  after  100 
cycles  were  obtained  using  CPFEM  and  illustrated  in  Figure  28. 


Figure  28:  Normalized  total  plastic  shear  strain  for  CP-Ti  (left)  and  Ti-0  (right). 
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As  clearly  seen  in  the  figure,  the  presence  of  oxygen  considerably  reduces  the  slip  activity  ahead 
of  the  notch  and  instead  facilitate  the  grain  boundary  decohesion  as  seen  previously.  Grain  level 
calculations  of  slips  can  be  used  to  estimate  the  probability  of  small  crack  initiation  and  growth 
within  the  microstructure. 

4.  Summary 

We  have  made  substantial  progress  towards  understanding  the  deformation  mechanism  by 
focusing  on  the  effect  of  grain  boundary  character  distributions  and  developing  the  necessary 
multiscale  tools.  Overall,  we  have  generalized  XFEM  to  include  implicit  representation  of  grains 
and  grain  boundaries.  The  new  formulation  can  incorporate  microstructure  features  directly  based 
on  EBSD  data  using  level  set  formulations.  Furthermore,  an  atomistic  database  of  interactions  of 
dislocation-impurities  was  developed.  Using  the  developed  database  an  appropriate  crystal 
plasticity  model  was  developed  to  investigate  SCG  in  Ti  alloys. 
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